//
// perc2 - calculation of the fractal dimension of correlated
// bond percolation cluster hulls
//
// Copyright (C) 2009, 2010 Indrek Mandre <indrek(at)mare.ee>
// http://www.mare.ee/indrek/perc2/, http://code.google.com/p/perc2/
// 
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
// 
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
// 
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
// 

/* Sample file:
A,M,SAMPLES,SUM(V1),SUM(V2),...,SUM(VM),SUM(V1*V1),SUM(V1*V2),...,SUM(V1*VM),SUM(V2*V2),SUM(V2*V3),...,SUM(VM*VM)
8,5,17179869184,430638803100,178192749907,25825988534,59757859615,42030918190,13387277869808,5282903952924,863072334900,1796212398546,1508124438734,2142236800867,326150657437,716937299364,554013470275,65462123698,106594146703,108534562097,339841259099,194128317063,203592808224
16,5,8589934592,743320637192,278631078459,47769220299,114480398760,95283581991,79551563705000,29049806433460,5203245160792,12158266070314,10899590155772,10719751140569,1881420588695,4455542107918,3881262302485,356529115457,776026911546,729649112234,2069137163218,1650212193731,1584758193491
32,5,4294967296,1264443890272,448569472136,82762634532,199897498996,184744733132,458954120745064,160506991554630,30177907059774,72421628237902,69275716577920,56358263001028,10527959125148,25353843309563,24012385941944,2011339074716,4727733298503,4579951221631,11807449637792,10905271416100,10664319328810
64,5,2147483648,2137247249996,735837807121,140624418511,340357734744,333325848286,2618662029633040,894739152628414,172517764919770,416819316279990,415118283508550,306155810609591,58906625921709,142461732216390,141402102459784,11410887284543,27401341695696,27386268325706,66999940250798,66031759133877,66233719753602
128,5,1073741824,3602726917480,1220944997144,237413808061,574952178005,580687411377,14873682757832864,5020617318339466,980506840047568,2373339575205436,2417119180378304,1695555652207068,330909166965631,801194805474931,815062220965243,64713617618243,156356279197657,159402876727379,379815382552775,385617833130098,393636729978811
256,5,536870912,6065721626954,2038932981578,399902520865,968616947460,994062127341,84307022680141068,28281270120357394,5558802963747540,13462221428058734,13873753993641190,9488679007027014,1864632493978563,4516101250862747,4652469441188009,366649527598227,887466525439064,914868515414199,2151534620149040,2215251032396837,2284651118251259
512,5,268435456,10206975447938,3416730285722,673019088779,1630216925238,1686824560275,477425626587527844,159650331154117826,31481070992807888,76251643680763644,79065031352869654,53389649341690290,10527053002232453,25498628013283292,26436373661341552,2076053502378859,5027690235932917,5213650466613613,12181628853905024,12627600975934367,13096559505791013
1024,5,134217728,17169616135924,5735362721153,1132160997309,2742416124706,2849478921599,2701779281043355432,902035054143748336,178156440647068622,431540215464416586,448858879233469126,301164860596834637,59480289594188940,144077400733140524,149854103894975424,11748073136185209,28455463291461809,29598213736036300,68932842180722864,71693392673719316,74576044898717029

Use it like this:
  $ ./calc -u 1 -u 2 -u 3 -s 3 test.csv

To visualize S(d) use the calc.gp:
  $ gnuplot calc.gp; epstopdf calc.eps; acroread calc.pdf
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <ctype.h>
#include <assert.h>

#ifdef __sun__
  #include <ieeefp.h>
#endif

#include <sstream>
#include <vector>
#include <deque>
#include <iostream>
#include <algorithm>

// Maximum rows of data allowed, increase if necessary
//#define MN 128
#define MN 9000
// maximum data columns allowed
#define MAXM 18

#define FTOL 3e-14

#if defined(PREC_LONGDOUBLE)
  typedef long double prec_t;
  static inline double v2d (const prec_t& v) { return v; } 
  static inline double mypow (const prec_t& x, const prec_t& y) { return powl(x, y); }
  static inline prec_t mylog (const prec_t& x) { return logl(x); }
#elif defined(PREC_GMP)
  #include <gmpxx.h>
  typedef mpf_class prec_t;
  static inline double v2d (const prec_t& v) { return v.get_d(); } 
  static inline prec_t mypow (const prec_t& x, const prec_t& y) { return pow(x.get_d(), y.get_d()); }
  static inline prec_t mylog (const prec_t& x) { return log(x.get_d()); }

  static std::istream& operator>> (std::istream& is, prec_t& obj)
  {
    double dbl;
    is >> dbl;
    obj = dbl;
    return is;
  }
#else
  typedef double prec_t;
  static inline double v2d (const prec_t& v) { return v; } 
  static inline double mypow (const prec_t& x, const prec_t& y) { return pow(x, y); }
  static inline prec_t mylog (const prec_t& x) { return log(x); }
  #undef PREC_DOUBLE
  #define PREC_DOUBLE
#endif

static inline prec_t myabs (const prec_t& x)
{
  if ( x < 0 )
      return -x;
  else
      return x;
}

static void error(const char *str)
{
  fprintf (stderr, "%s\n", str);
  throw str;
}

class matrix
{
public:
  matrix() : rows(0), cols(0) { }

  matrix(size_t rows, size_t cols) : rows(rows), cols(cols)
  {
    if ( !rows || !cols ) throw "can't have nonzero count of rows or columns in a matrix";
    cells.resize(rows*cols);
  }

  matrix(size_t rows, size_t cols, prec_t value) : rows(rows), cols(cols)
  {
    if ( !rows || !cols ) throw "can't have nonzero count of rows or columns in a matrix";
    cells.resize(rows*cols, value);
  }

  matrix(size_t rows, size_t cols, const prec_t *vals) : rows(rows), cols(cols)
  {
    if ( !rows || !cols ) throw "can't have nonzero count of rows or columns in a matrix";
    cells.resize(rows*cols);
    for ( size_t r = 0; r < rows; r++ )
        for ( size_t c = 0; c < cols; c++ )
            at(r, c) = *vals++;
  }

  matrix(const matrix& m) : rows(m.rows), cols(m.cols), cells(m.cells) { }

  void clear ()
  {
    cells.clear();
    rows = cols = 0;
  }

  inline const prec_t& at(size_t i, size_t j) const { return cells[i * cols + j]; }
  inline prec_t& at(size_t i, size_t j) { return cells[i * cols + j]; }

  inline const prec_t& operator() (size_t i, size_t j) const { return at(i, j); }
  inline prec_t& operator() (size_t i, size_t j) { return at(i, j); }

  inline void mul_row(size_t row, const prec_t& v)
  {
    for ( size_t i = 0; i < cols; i++ )
        (*this)(row, i) *= v;
  }

  inline void mul_col(size_t col, const prec_t& v)
  {
    for ( size_t i = 0; i < rows; i++ )
        (*this)(i, col) *= v;
  }

  inline void add_row(size_t row, size_t from, const prec_t& mul)
  {
    for ( size_t j = 0; j < cols; j++ )
        at(row, j) += mul * at(from, j);
  }

  inline void add_col(size_t col, size_t from, const prec_t& mul)
  {
    for ( size_t j = 0; j < rows; j++ )
        (*this)(j, col) += mul * (*this)(j, from);
  }

  inline void swap_rows(size_t r1, size_t r2)
  {
    for ( size_t i = 0; i < cols; i++ )
        std::swap((*this)(r1, i), (*this)(r2, i));
  }

  inline void swap_cols(size_t c1, size_t c2)
  {
    for ( size_t i = 0; i < rows; i++ )
        std::swap((*this)(i, c1), (*this)(i, c2));
  }

  const matrix& operator= (const matrix& m)
  {
    cols = m.cols;
    rows = m.rows;
    cells = m.cells;
    return *this;
  }

  prec_t det() const
  {
    throw("todo");
    return 0;
  }

  void invert ()
  {
    if ( rows != cols || !rows ) throw "we need a nonzero square matrix to pivot!";
    matrix m(cols, rows, 0.0);
    for ( size_t i = 0; i < cols; i++ )
        m(i, i) = 1.0;
    for ( size_t c = 0; c < cols; c++ )
      {
        // first find largest pivot in this column
        size_t mpr = c;
        for ( size_t r = c + 1; r < rows; r++ )
            if ( myabs(at(mpr, c)) < myabs(at(r, c)) )
                mpr = r;
        if ( mpr != c )
          {
            swap_rows(mpr, c);
            m.swap_rows(mpr, c);
          }
        prec_t pivot = at(c, c);
        // now reduce rows above and below this to zero
        for ( size_t r = 0; r < rows; r++ )
          {
            if ( r == c )
                continue;
            prec_t f =  -at(r, c) / pivot;
            add_row(r, c, f);
            m.add_row(r, c, f);
          }
        // now make this row into 1-s
        mul_row(c, 1.0 / pivot);
        m.mul_row(c, 1.0 / pivot);
      }
    // we're done, swap the inverted matrix from the right side
    cells.swap(m.cells);
  }

  void print () const
  {
    for ( size_t r = 0; r < rows; r++ )
      {
        for ( size_t c = 0; c < cols; c++ )
          {
            if ( at(r,c) >= 0 )
                printf (" ");
            printf (" %.5f", v2d(at(r, c)));
          }
        printf ("\n");
      }
  }

  static void mul(matrix& m, const matrix& lm, const matrix& rm)
  {
    if ( lm.cols != rm.rows ) throw "incompatible sized matrices for multiplication";
    m.resize(lm.rows, rm.cols);
    for ( size_t r = 0; r < m.rows; r++ )
      {
        for ( size_t c = 0; c < m.cols; c++ )
          {
            m(r, c) = 0;
            for ( size_t i = 0; i < lm.cols; i++ )
                m(r, c) += lm(r, i) * rm(i, c);
          }
      }
  }

  void resize(size_t r, size_t c)
  {
    rows = r;
    cols = c;
    if ( !rows || !cols ) throw "can't have nonzero count of rows or columns in a matrix";
    cells.resize(r * c);
  }

  void resize(size_t r, size_t c, const prec_t& value)
  {
    rows = r;
    cols = c;
    if ( !rows || !cols ) throw "can't have nonzero count of rows or columns in a matrix";
    cells.resize(r * c, value);
  }

  void swap(matrix& m)
  {
    std::swap(rows, m.rows);
    std::swap(cols, m.cols);
    cells.swap(m.cells);
  }

private:
  size_t rows, cols;
  std::vector<prec_t> cells;
};

// Generated using Octave with:
//  b = cat(2, (1:128)', chi2inv (0.95, (1:128)'));
//  printf ("  ");
//  for i=1:size(b, 1)
//    printf ("{%3d, %12.8f},", b(i,1), b(i,2));
//    if mod(i,4) == 0
//      printf ("\n  ");
//    else
//      printf (" ");
//    endif
//  endfor

static const struct {
  int dof;
  prec_t val;
} chi2inv_data[] = {
  {  0,   0         },
  {  1,   3.84145882}, {  2,   5.99146455}, {  3,   7.81472790}, {  4,   9.48772904},
  {  5,  11.07049769}, {  6,  12.59158724}, {  7,  14.06714045}, {  8,  15.50731306},
  {  9,  16.91897760}, { 10,  18.30703805}, { 11,  19.67513757}, { 12,  21.02606982},
  { 13,  22.36203249}, { 14,  23.68479130}, { 15,  24.99579014}, { 16,  26.29622760},
  { 17,  27.58711164}, { 18,  28.86929943}, { 19,  30.14352721}, { 20,  31.41043284},
  { 21,  32.67057334}, { 22,  33.92443847}, { 23,  35.17246163}, { 24,  36.41502850},
  { 25,  37.65248413}, { 26,  38.88513866}, { 27,  40.11327207}, { 28,  41.33713815},
  { 29,  42.55696780}, { 30,  43.77297183}, { 31,  44.98534328}, { 32,  46.19425952},
  { 33,  47.39988392}, { 34,  48.60236737}, { 35,  49.80184957}, { 36,  50.99846017},
  { 37,  52.19231973}, { 38,  53.38354062}, { 39,  54.57222776}, { 40,  55.75847928},
  { 41,  56.94238715}, { 42,  58.12403768}, { 43,  59.30351203}, { 44,  60.48088658},
  { 45,  61.65623338}, { 46,  62.82962041}, { 47,  64.00111197}, { 48,  65.17076890},
  { 49,  66.33864886}, { 50,  67.50480655}, { 51,  68.66929391}, { 52,  69.83216034},
  { 53,  70.99345283}, { 54,  72.15321617}, { 55,  73.31149303}, { 56,  74.46832416},
  { 57,  75.62374847}, { 58,  76.77780316}, { 59,  77.93052381}, { 60,  79.08194449},
  { 61,  80.23209785}, { 62,  81.38101519}, { 63,  82.52872654}, { 64,  83.67526074},
  { 65,  84.82064550}, { 66,  85.96490744}, { 67,  87.10807220}, { 68,  88.25016442},
  { 69,  89.39120787}, { 70,  90.53122543}, { 71,  91.67023918}, { 72,  92.80827038},
  { 73,  93.94533960}, { 74,  95.08146667}, { 75,  96.21667075}, { 76,  97.35097038},
  { 77,  98.48438346}, { 78,  99.61692732}, { 79, 100.74861875}, { 80, 101.87947397},
  { 81, 103.00950871}, { 82, 104.13873823}, { 83, 105.26717730}, { 84, 106.39484024},
  { 85, 107.52174097}, { 86, 108.64789297}, { 87, 109.77330935}, { 88, 110.89800282},
  { 89, 112.02198575}, { 90, 113.14527014}, { 91, 114.26786768}, { 92, 115.38978971},
  { 93, 116.51104728}, { 94, 117.63165114}, { 95, 118.75161175}, { 96, 119.87093930},
  { 97, 120.98964370}, { 98, 122.10773461}, { 99, 123.22522145}, {100, 124.34211340},
  {101, 125.45841941}, {102, 126.57414819}, {103, 127.68930826}, {104, 128.80390793},
  {105, 129.91795529}, {106, 131.03145825}, {107, 132.14442454}, {108, 133.25686170},
  {109, 134.36877710}, {110, 135.48017793}, {111, 136.59107123}, {112, 137.70146386},
  {113, 138.81136256}, {114, 139.92077390}, {115, 141.02970429}, {116, 142.13816004},
  {117, 143.24614728}, {118, 144.35367205}, {119, 145.46074022}, {120, 146.56735758},
  {121, 147.67352976}, {122, 148.77926230}, {123, 149.88456062}, {124, 150.98943002},
  {125, 152.09387569}, {126, 153.19790274}, {127, 154.30151617}, {128, 155.40472085}
};

static const prec_t& chi2inv (size_t n)
{
  if ( n > 128 )
    {
      fprintf (stderr, "chi2inv: table overflow with %zu\n", n);
      static prec_t zero = 0;
      return zero;
      exit(-1);
    }
  return chi2inv_data[n].val;
}

static size_t N = 0;
static size_t M = 0;
static prec_t A[MN];
static prec_t NS[MN];
static prec_t X[MN][MAXM];
static prec_t COVM[MN][MAXM][MAXM];
static prec_t ICOVM[MN][MAXM][MAXM];
static bool EXTRA_MIN = false;

void calc_icovm()
{
  matrix m(M, M);
  for ( size_t i = 0; i < N; i++ )
    {
      for ( size_t j = 0; j < M; j++ )
          for ( size_t k = 0; k < M; k++ )
              m(j,k) = COVM[i][j][k];
      m.invert();
      for ( size_t j = 0; j < M; j++ )
          for ( size_t k = 0; k < M; k++ )
              ICOVM[i][j][k] = m(j,k);
    }
}

struct chicalc
{
public:
  inline chicalc(const prec_t *Y) : Y(Y) { }

  prec_t operator() (const std::vector<prec_t>& P)
  {
    prec_t S = 0;
    for ( size_t i = 0; i < N; i++ )
      {
        prec_t m = 0;
        for ( size_t k = 0; k < M; k++ )
          m += P[k] * X[i][k];
        prec_t s2 = 0;
        for ( size_t k = 0; k < M; k++ )
            for ( size_t l = 0; l < M; l++ )
                s2 += P[k]*P[l]*COVM[i][k][l];
        S += (Y[i] - m) * (Y[i] - m) / s2;
      }
    return S;
  }

  void df(const std::vector<prec_t>& P, std::vector<prec_t>& DF)
  {
    prec_t top[N];
    prec_t top2[N];
    prec_t bottom[N];
    prec_t bottom2[N];
    for ( size_t i = 0; i < N; i++ )
      {
        prec_t m = 0;
        for ( size_t k = 0; k < M; k++ )
          m += P[k] * X[i][k];
        bottom[i] = 0;
        for ( size_t k = 0; k < M; k++ )
            for ( size_t l = 0; l < M; l++ )
                bottom[i] += P[k]*P[l]*COVM[i][k][l];
        bottom2[i] = bottom[i] * bottom[i];
        top[i] = Y[i] - m;
        top2[i] = top[i] * top[i];
      }
    for ( size_t c = 0; c < M; c++ )
      {
        prec_t dS = 0;
        for ( size_t i = 0; i < N; i++ )
          {
            prec_t d_top = -2 * top[i] * X[i][c];
            prec_t d_bottom = 0;
            for ( size_t k = 0; k < M; k++ )
              {
                for ( size_t l = 0; l < M; l++ )
                    {
                      if ( k == c && l == c )
                          d_bottom += 2 * P[l]*COVM[i][k][l];
                      else if ( k == c )
                          d_bottom += P[l]*COVM[i][k][l];
                      else if ( l == c )
                          d_bottom += P[k]*COVM[i][k][l];
                    }
              }
            dS += (d_top * bottom[i] - d_bottom * top2[i]) / bottom2[i];
          }
        DF[c] = dS;
      }
  }

private:
  const prec_t *Y;
};

static std::vector<prec_t> C;

static inline prec_t SIGN(const prec_t& l, const prec_t& r)
{
  if ( r >= 0 )
      return myabs(l);
  else
      return -myabs(r);
}

static inline void shft3(prec_t &a, prec_t& b, prec_t& c, const prec_t& d)
{
  a = b;
  b = c;
  c = d;
}

struct bracket
{
  prec_t fa, fb, fc;
  prec_t ax, bx, cx;

  bracket() { }

  template <class F>
  void run(F& func, const prec_t& a, const prec_t& b)
  {
    const prec_t GS = 1.618034;
    const prec_t GSLIMIT = 100.0;
    const prec_t TINY = 1.0e-20;
    ax = a;
    bx = b;
    fa = func(ax);
    fb = func(bx);
    if ( fb > fa )
      {
        std::swap(ax, bx);
        std::swap(fa, fb);
      }
    cx = bx + GS * (bx - ax);
    fc = func(cx);

    while ( fb > fc )
      {
        prec_t r = (bx - ax) * (fb - fc);
        prec_t q = (bx - cx) * (fb - fa);
        prec_t u = bx - ((bx - cx) * q - (bx - ax) * r) / (2.0 *  SIGN(std::max(myabs(q - r), TINY), q - r));
        prec_t ulim = bx + GSLIMIT * (cx - bx);
        prec_t fu;
        if ( (bx - u) * (u - cx) > 0.0 )
          {
            fu = func(u);
            if ( fu < fc )
              {
                ax = bx;
                bx = u;
                fa = fb;
                fb = fu;
                return;
              }
            else if ( fu > fb )
              {
                cx = u;
                fc = fu;
                return;
              }
            u = cx + GS * (cx - bx);
            fu = func(u);
          }
        else if ( (cx - u) * (u - ulim) > 0.0 )
          {
            fu = func(u);
            if ( fu < fc )
              {
                shft3(bx, cx, u, u + GS * (u - cx));
                shft3(fb, fc, fu, func(u));
              }
          }
        else if ( (u - ulim) * (ulim - cx) >= 0.0 )
          {
            u = ulim;
            fu = func(u);
          }
        else
          {
            u = cx + GS * (cx - bx);
            fu = func(u);
          }
        shft3(ax, bx, cx, u);
        shft3(fa, fb, fc, fu);
      }
  }

};

template <class F>
struct F1dim
{
  F func;
  typedef std::vector<prec_t> V;
  V start, dir;
  inline F1dim(const F& func, const V& start, const V& dir) : func(func), start(start), dir(dir) { }

  prec_t operator() (const prec_t& t)
  {
    V p = start;
    for ( size_t i = 0; i < p.size(); i++ )
        p[i] += t * dir[i];
    return func(p);
  }
};

struct Brent
{
public:
  prec_t xmin, fmin;
  const prec_t tol;
  Brent(const prec_t& tol = FTOL) : tol(tol) { }

  template <class T>
  prec_t minimize (T& func, const prec_t& ax, const prec_t& bx, const prec_t& cx)
  {
    const size_t ITMAX = 2000;
    const prec_t CGOLD = 0.3819660;
    const prec_t ZEPS = 3e-21;

    prec_t a, b, d = 0, etemp, fu, fv, fw, fx;
    prec_t p, q, r, tol1, tol2, u, v, w, x, xm;
    prec_t e = 0.0;
    a = ax;
    b = cx;
    if ( a > b )
        std::swap(a, b);
    x = w = v = bx;
    fw = fv = fx = func(x);
    for ( size_t iter = 0; iter < ITMAX; iter++ )
      {
        xm = 0.5 * (a + b);
        tol2 = 2.0 * (tol1 = tol * myabs(x) + ZEPS);
        if ( myabs(x - xm) <= (tol2 - 0.5 * (b - a)))
          {
            fmin = fx;
            xmin = x;
            return xmin;
          }
        if ( myabs(e) > tol1 )
          {
            r = (x - w) * (fx - fv);
            q = (x - v) * (fx - fw);
            p = (x - v) * q - (x - w ) * r;
            q = 2.0 * (q - r);
            if ( q > 0.0 )
                p = -p;
            q = myabs(q);
            etemp = e;
            e = d;
            if ( myabs(p) >= myabs ( 0.5 * q * etemp) || p <= q * (a - x) || p >= q * (b - x) )
              {
                d = CGOLD * (e = (x >= xm ? a - x : b - x));
              }
            else
              {
                d = p / q;
                u = x + d;
                if ( u - a < tol2 || b - u < tol2 )
                    d = SIGN(tol1, xm - x);
              }
          }
        else
          {
            d = CGOLD * (e=(x >= xm ? a - x : b - x));
          }
        u = (myabs(d) >= tol1 ? x + d : x + SIGN(tol1, d));
        fu = func(u);
        if ( fu <= fx )
          {
            if ( u >= x )
                a = x;
            else 
                b = x;
            shft3(v, w, x, u);
            shft3(fv, fw, fx, fu);
          }
        else
          {
            if ( u < x )
                a = u;
            else
                b = u;
            if ( fu <= fw || w == x )
              {
                v = w;
                w = u;
                fv = fw;
                fw = fu;
              }
            else if ( fu <= fv || v == x || v == w )
              {
                v = u;
                fv = fu;
              }
          }
      }
    fprintf (stderr, "too many iterations in brent!\n");
    fmin = fx;
    xmin = x;
    return xmin;
  }
};

template <class T>
struct Linemethod
{
  std::vector<prec_t> p;
  std::vector<prec_t> xi;
  T& func;
  size_t n;
  Linemethod(T& func) : func(func) { }

  prec_t linmin()
  {
    prec_t xmin;
    n = p.size();

    F1dim<T> f1dim(func, p, xi);
    bracket bra;
    bra.run(f1dim, 0.0, 1);
    Brent brent;
    xmin = brent.minimize(f1dim, bra.ax, bra.bx, bra.cx);
    for ( size_t j = 0; j < n; j++ )
      {
        xi[j] *= xmin;
        p[j] += xi[j];
      }
    return brent.fmin;
  }
};

template <class T>
struct Frprmn : Linemethod<T>
{
  size_t iter;
  prec_t fret;
  using Linemethod<T>::func;
  using Linemethod<T>::linmin;
  using Linemethod<T>::p;
  using Linemethod<T>::xi;
  const prec_t ftol;

  Frprmn(T& func, const prec_t& ftol = FTOL) : Linemethod<T>(func), ftol(ftol) { }

  prec_t minimize(const std::vector<prec_t>& pp)
  {
    const size_t ITMAX = 2000;
    const prec_t EPS = 1.0e-18;
    const prec_t GTOL = 1.0e-8;
    size_t countdown = 3;
    prec_t gg, dgg;
    size_t n = pp.size();
    p = pp;
    std::vector<prec_t> g(n), h(n);
    xi.resize(n);
    prec_t fp = func(p);
    func.df(p, xi);
    for ( size_t j = 0; j < n; j++ )
      {
        g[j] = -xi[j];
        xi[j] = h[j] = g[j];
      }
    for ( size_t its = 0; its < ITMAX; its++ )
      {
        iter = its;
        std::vector<prec_t> tmp = p;
        fret = linmin();
        if ( 2.0 * myabs(fret - fp) <= ftol * (myabs(fret) + myabs(fp) + EPS) && !--countdown )
            return fret;
        fp = fret;
        func.df(p, xi);
        prec_t test = 0.0;
        prec_t den = std::max(fp, (prec_t)1.0);
        for ( size_t j = 0; j < n; j++ )
          {
            prec_t temp = myabs(xi[j]) * std::max(myabs(p[j]), (prec_t)1.0) / den;
            if ( temp > test )
                test = temp;
          }
        if ( test < GTOL )
            return fret;
        dgg = gg = 0.0;
        for ( size_t j = 0; j < n; j++ )
          {
            gg += g[j] * g[j];
            dgg += (xi[j] + g[j]) * xi[j];
          }
        if ( gg == 0.0 )
            return fret;
        prec_t gam = dgg / gg;
        for ( size_t j = 0; j < n; j++ )
          {
            g[j] = -xi[j];
            xi[j] = h[j] = g[j] + gam * h[j];
          }
      }
    //error ("Too many iterations in frprmn!");
    return fret;
  }
};


prec_t optimize(const std::vector<prec_t>& P, chicalc& cc, std::vector<prec_t>& C, prec_t& oldS)
{
  Frprmn<chicalc> frp(cc);
  prec_t newS = frp.minimize(P);
  if ( newS < oldS )
    {
      C = frp.p;
      oldS = newS;
    }
  return newS;
}

// Calculate the residuals at the given exponent
static prec_t calc_res (prec_t d, prec_t *retC = 0)
{
  static bool first = true;
  if ( first )
    {
      C.resize(M, 1.0);
      first = false;
    }

  prec_t Y[N];
  for ( size_t i = 0; i < N; i++ )
      Y[i] = mypow(A[i], d);

  chicalc cc(Y);
  prec_t S = 9e33;

  // find result using classical least squares (iteratively)
  std::vector<prec_t> LC(M, 1.0);
  for ( size_t ci = 0; ci < 10; ci++ )
    {
      prec_t W[N];
      for ( size_t i = 0; i < N; i++ )
        {
          W[i] = 0;
          for ( size_t k = 0; k < M; k++ )
              for ( size_t l = 0; l < M; l++ )
                  W[i] += LC[k]*LC[l]*COVM[i][k][l];
          W[i] = 1/W[i];
        }

      prec_t XY[M][M + 1];
      for ( size_t k = 0; k < M; k++ )
        {
          for ( size_t l = 0; l < M; l++ )
            {
              XY[k][l] = 0;
              for ( size_t i = 0; i < N; i++ )
                  XY[k][l] += W[i] * X[i][k] * X[i][l];
            }
        }

      prec_t YW[N];
      for ( size_t i = 0; i < N; i++ )
          YW[i] = W[i] * Y[i];

      for ( size_t j = 0; j < M; j++ )
        {
          XY[j][M] = 0;
          for ( size_t i = 0; i < N; i++ )
              XY[j][M] += X[i][j] * YW[i];
        }

      // pivoting array
      prec_t *P[M];
      for ( size_t i = 0; i < M; i++ )
          P[i] = XY[i];

      for ( size_t i = 0; i < M; i++ )
        {
          // find largest pivot and move it to front
          size_t maxp = i;
          for ( size_t j = i + 1; j < M; j++ )
              if ( myabs(P[maxp][i]) < myabs(P[j][i]) )
                  maxp = j;
          std::swap (P[i], P[maxp]);

          // eliminate
          for ( size_t j = i + 1; j < M; j++ )
            {
              prec_t mod = P[j][i] / P[i][i];
              for ( size_t k = i; k <= M; k++ )
                  P[j][k] = P[j][k] - mod * P[i][k];
            }
        }

      // backsubstitute and get the results
      for ( size_t i = 0; i < M; i++ )
        {
          size_t idx = M - i - 1;
          LC[idx] = P[idx][M];
          for ( size_t j = idx + 1; j < M; j++ )
              LC[idx] -= P[idx][j] * LC[j];
          LC[idx] = LC[idx] / P[idx][idx];
        }
    }

  if ( EXTRA_MIN )
    {
      // first try to move from olde input
      optimize(C, cc, C, S);

      // then try from traditional input
      optimize(LC, cc, C, S);

      // now probe directionally along gully of minimum
      std::vector<prec_t> dir(M);
      // find offset amount
      prec_t mx = 0;
      for ( size_t i = 0; i < M; i++ )
          mx = std::max(mx, myabs(C[i]));
      // follow gradient to bottom
      Linemethod<chicalc> lm(cc);
      lm.p.resize(M);
      for ( size_t i = 0; i < M; i++ )
          lm.p[i] = C[i];
      lm.xi.resize(M);
      // offset current position
      lm.p[0] += std::max(mx / 100.0, (prec_t)1.0);
      cc.df(lm.p, lm.xi);
      lm.linmin();
      // now get a direction, normalize it
      prec_t l2 = 0;
      for ( size_t i = 0; i < M; i++ )
        {
          dir[i] = lm.p[i] - C[i];
          l2 += dir[i] * dir[i];
        }
      prec_t l = sqrt(l2);
      for ( size_t i = 0; i < M; i++ )
          dir[i] = dir[i] / l;

      // now try to optimize in either direction
      std::vector<prec_t> w(M);
      for ( size_t i = 0; i < M; i++ )
          w[i] = LC[i] + dir[i] * 1;
      optimize(w, cc, C, S);
      for ( size_t i = 0; i < M; i++ )
          w[i] = LC[i] + -dir[i] * 1;
      optimize(w, cc, C, S);
    }
  else
    {
      C = LC;
      S = cc(C);
    }

  if ( !finite (v2d(S)) )
    {
      printf ("NaN/Inf at d=%.8f, returning 1e4 instead\n", v2d(d));
      return 1e4;
    }
  if ( retC )
    {
      for ( size_t i = 0; i < M; i++ )
          retC[i] = C[i];
    }

#if 0
  prec_t mm = 1e99;
  FILE *fp = fopen("map.csv", "w");
  #define MRANGE 2000
  printf ("old=%g %g\n", v2d(C[0]), v2d(C[1]));
  const prec_t add = 1000;
  prec_t c0s = C[0] - add;
  prec_t c0e = C[0] + add;
  prec_t c1s = C[1] - add;
  prec_t c1e = C[1] + add;
  for ( size_t j = 0; j < MRANGE; j++ )
    {
      C[1] =  c1s + j * (c1e - c1s) / (MRANGE - 1);
      for ( size_t i = 0; i < MRANGE; i++ )
        {
          C[0] =  c0s + i * (c0e - c0s) / (MRANGE - 1);
          prec_t S = calc_chi(Y, C);
          fprintf (fp, "%g ", v2d(S));
          if ( mm > S )
            {
              //printf ("best=%g %g\n", v2d(C[0]), v2d(C[1]));
              mm = S;
            }
        }
      fprintf (fp, "\n");
    }
  fclose (fp);
  printf ("min=%g\n", v2d(mm));
  exit(0);
#endif

  return S;
}

struct calc_res_functor
{
  inline prec_t operator()(const prec_t& v) const { return calc_res(v); }
};

static void showC (prec_t d)
{
  prec_t C[M];
  calc_res(d, C);

  printf ("A^%g = ", v2d(d));
  for ( size_t i = 0; i < M; i++ )
    {

      if ( i != 0 ) printf (" + ");
      printf ("%.8e * X%zu", v2d(C[i]), i + 1);
    }
  printf ("\n");
}

// Find a minimum using the golden search method
template <class Functor>
static prec_t find_min (Functor F, prec_t a, prec_t b, prec_t c)
{
  prec_t ra = F (a);
  prec_t rb = F (b);
  prec_t rc = F (c);

#ifdef PREC_DOUBLE
  while ( c - a > 1e-8 )
    {
      if ( (ra > rc && b - a > 1e-8 / 2) || c - b < 1e-8 / 2 )
#else
  while ( c - a > 1e-16 )
    {
      if ( (ra > rc && b - a > 1e-16 / 2) || c - b < 1e-16 / 2 )
#endif
        {
          prec_t m = 0.5 * (a + b);
          prec_t rm = F (m);
          if ( rm > rb )
            {
              a = m;
              ra = rm;
            }
          else
            {
              c = b;
              rc = rc;
              b = m;
              rb = rm;
            }
        }
      else
        {
          prec_t m = 0.5 * (b + c);
          prec_t rm = F (m);
          if ( rm > rb )
            {
              c = m;
              rc = rm;
            }
          else
            {
              a = b;
              ra = rb;
              b = m;
              rb = rm;
            }
        }
    }

  return 0.5 * (a + c);
}

struct find_value_min
{
  find_value_min (const prec_t& sub) : sub(sub) { }

  inline prec_t operator() (const prec_t& x)
  {
    prec_t ret = calc_res (x) - sub;
    if ( ret < 0 )
        return -ret;
    else
        return ret;
  }

  prec_t sub;
};

#define FIND_VALUE_DIVF 10

static prec_t find_value (prec_t value, prec_t start, prec_t dir, prec_t limit)
{
  find_value_min fv(value);
  prec_t add = dir / FIND_VALUE_DIVF;
  bool posit = dir > 0;

  // adjust our picked increment accordingly
  for ( size_t i = 0; i < 10; i++ )
    {
      // do this until the next point is better than the start
      prec_t v1 = fv(start);
      prec_t v2 = fv(start + add);
      prec_t v3 = fv(start + 3 * add);
      if ( v1 > v2 && v2 > v3 )
          break;
      if ( i == 9 ) // data does not make sense
          return start;
      add = add / 2;
    }

  prec_t x = start;
  prec_t minval = fv(x);
  prec_t minx = x;

  // step on and find the area of the minimum
  while ( 1 )
    {
      x += add;
      if ( (posit && x > limit) || (!posit && x < limit) ) // failure, out of bounds
          return start;
      prec_t y = fv(x);
      if ( y < minval )
        {
          minval = y;
          minx = x;
        }
      else
        {
          break;
        }
      add *= 1.25;
      if ( (posit && add > dir) || (!posit && add < dir) )
          add = dir;
    }

  // use golden search to find the minimum
  if ( posit )
      return find_min (fv, start, minx, x);
  else
      return find_min (fv, x, minx, start);
}

static void usage (const char *argv0)
{
  fprintf (stderr, "%s: %s [-r low high] [-i intervals] [-s skip] [-d discard [-d discard ..]] [-f exponent [-f exponent ..]] [-a column exponent amplitude] [-p] [data.csv]\n", argv0, argv0);
  fprintf (stderr, "    -r low high - specify the search range\n");
  fprintf (stderr, "    -i intervals - divide the search range into given number of intervals\n");
  fprintf (stderr, "    -s skip - skip given number of data records\n");
  fprintf (stderr, "    -d discard - discard given data column\n");
  fprintf (stderr, "    -u use - use given data column\n");
  fprintf (stderr, "    -f exponent - add a fixed exponent\n");
  fprintf (stderr, "    -a column exponent amplitude - adds a static component to the specified column\n");
  fprintf (stderr, "    -p - display processed data and covariance matrices\n");
  fprintf (stderr, "    -t - apply traditional chi2 minimization for all data following the algorithm\n");
  fprintf (stderr, "    -l - further minimize the transformed chi2 using gradient\n");
  fprintf (stderr, "    -n - do not calculate minimums\n");
  exit (-1);
}

static matrix RM;
static std::vector<prec_t> D;

static matrix alpha, beta;

prec_t calc_y(size_t n, size_t m)
{
  prec_t res = 0;
  for ( size_t r = 0; r < M; r++ )
        res += RM(m, r)*mypow(A[n], D[r]);
  return res;
}

prec_t calc_dydd(size_t n, size_t m, size_t di)
{
  return RM(m, di)*mylog(A[n])*mypow(A[n], D[di]);
}

prec_t calc_dydf(size_t n, size_t m, size_t fr, size_t fc)
{
  if ( fr != m )
      return 0;
  return mypow(A[n], D[fc]);
}

struct x_c
{
  x_c() { }
  x_c(size_t n, size_t m) : n(n), m(m) { }

  size_t n;
  size_t m;
};

#define INVALID_VALUE ((size_t)-1)

struct a_idx
{
  a_idx() { }
  a_idx(size_t di, size_t r, size_t c) : di(di), r(r), c(c) { }
  size_t di;
  size_t r, c;
};

prec_t calc_y(const x_c& coord)
{
  prec_t res = 0;
  for ( size_t r = 0; r < M; r++ )
        res += RM(coord.m, r)*mypow(A[coord.n], D[r]);
  return res;
}

prec_t calc_dyda(const x_c& coord, struct a_idx& aid)
{
  if ( aid.di != INVALID_VALUE ) // derivative against exponent
    {
      return RM(coord.m, aid.di)*mylog(A[coord.n])*mypow(A[coord.n], D[aid.di]);
    }
  else
    {
      if ( aid.r != coord.m )
          return 0;
      return mypow(A[coord.n], D[aid.c]);
    }
}

void make_alpha_beta ()
{
  std::vector<prec_t> y;
  std::vector<x_c> coords;
  for ( size_t i = 0; i < N; i++ )
    {
      for ( size_t j = 0; j < M; j++ )
        {
          coords.push_back (x_c(i, j));
          y.push_back (X[i][j]);
        }
    }

  matrix W(coords.size(), coords.size(), 0.0);
  for ( size_t i = 0; i < coords.size(); i++ )
    {
      for ( size_t j = 0; j < coords.size(); j++ )
        {
          if ( coords[i].n != coords[j].n )
              continue;
          W(i, j) = ICOVM[coords[i].n][coords[i].m][coords[j].m];
        }
    }

  std::vector<a_idx> params;
  for ( size_t j = 0; j < M; j++ )
      params.push_back (a_idx(j, INVALID_VALUE, INVALID_VALUE));
  for ( size_t j = 0; j < M; j++ )
      for ( size_t k = 0; k < M; k++ )
          params.push_back (a_idx(INVALID_VALUE, j, k));

  beta.clear();
  beta.resize(params.size(), 1, 0);
  for ( size_t k = 0; k < params.size(); k++ )
    {
      beta(k, 0) = 0.0;
      for ( size_t i = 0; i < coords.size(); i++ )
        {
          for ( size_t j = 0; j < coords.size(); j++ )
            {
              beta(k, 0) += 0.5 * W(i, j) * (calc_dyda(coords[i], params[k]) * (y[j] - calc_y(coords[j])) +
                  calc_dyda(coords[j], params[k]) * (y[i] - calc_y(coords[i])));
                  
            }
        }
    }
  //printf ("Beta:\n");
  //beta.print();

  alpha.clear();
  alpha.resize(params.size(), params.size(), 0);
  for ( size_t k = 0; k < params.size(); k++ )
    {
      for ( size_t l = 0; l < params.size(); l++ )
        {
          alpha(k, l) = 0.0;
          for ( size_t i = 0; i < coords.size(); i++ )
            {
              for ( size_t j = 0; j < coords.size(); j++ )
                {
                  alpha(k, l) += 0.5 * W(i, j) * (calc_dyda(coords[i], params[k]) * calc_dyda(coords[j], params[l]) + calc_dyda(coords[j], params[k]) * calc_dyda(coords[i], params[l]));
                      
                }
            }
        }
    }
  //alpha.print();
}

prec_t calc_chi2overall()
{
#if 0
  prec_t sum = 0.0;
  for ( size_t i = 0; i < N; i++ )
    {
      for ( size_t j = 0; j < M; j++ )
        {
          prec_t res = calc_y(i, j);
          prec_t diff = res - X[i][j];
          //printf ("%g-%g=%g, stdev=%g\n", v2d(res), v2d(X[i][j]), v2d(diff), v2d(sqrt(COVM[i][j][j])));
          sum += diff * diff / COVM[i][j][j];
        }
    }
#endif
  prec_t sum = 0.0;
  for ( size_t i = 0; i < N; i++ )
    {
      for ( size_t j = 0; j < M; j++ )
        {
          for ( size_t k = 0; k < M; k++ )
            {
              prec_t w_jk = ICOVM[i][j][k];
              prec_t res1 = calc_y(i, j);
              prec_t diff1 = res1 - X[i][j];
              prec_t res2 = calc_y(i, k);
              prec_t diff2 = res2 - X[i][k];
              sum += diff1 * diff2 * w_jk;
            }
          //printf ("%g-%g=%g, stdev=%g\n", v2d(res), v2d(X[i][j]), v2d(diff), v2d(sqrt(COVM[i][j][j])));
        }
    }
  return sum;
}

int main (int argc, char *argv[])
{
  double low = -0.5;
  double high = 2.0;
  size_t intervals = 5000;
  size_t skip = 0;
  std::vector<size_t> dv;
  std::vector<size_t> uv;
  std::vector<double> fv;
  size_t add_column = 0;
  double add_exponent;
  double add_amplitude;
  bool disp = false;
  bool apply_traditional = false;
  bool nomins = false;

  const char *fn = "calc.csv";

  for ( int i = 1; i < argc; i++ )
    {
      if ( argv[i][0] != '-' )
        {
          fn = argv[i];
          continue;
        }

      if ( argv[i][0] == '-' && (argv[i][1] == 0 || argv[i][2] != 0) )
          usage (argv[0]);

      switch (argv[i][1])
        {
          case 'n': {
            nomins = true;
            break;
          };
          case 't': {
            apply_traditional = true;
            break;
          };
          case 'p': {
            disp = true;
            break;
          };
          case 'f': {
            double fixed;
            if ( i + 1 >= argc || sscanf (argv[++i], "%lf", &fixed) != 1 )
                usage (argv[0]);
            if ( std::find(fv.begin(), fv.end(), fixed) != fv.end() )
                usage (argv[0]);
            fv.push_back (fixed);
            break;
          };
          case 'd': {
            size_t discard;
            if ( i + 1 >= argc || sscanf (argv[++i], "%zu", &discard) != 1 || !uv.empty() )
                usage (argv[0]);
            if ( std::find(dv.begin(), dv.end(), discard) != dv.end() )
                usage (argv[0]);
            dv.push_back (discard);
            break;
          };
          case 'u': {
            size_t use;
            if ( i + 1 >= argc || sscanf (argv[++i], "%zu", &use) != 1 || !dv.empty() )
                usage (argv[0]);
            if ( std::find(uv.begin(), uv.end(), use) != uv.end() )
                usage (argv[0]);
            uv.push_back (use);
            break;
          };
          case 's': {
            if ( i + 1 >= argc || sscanf (argv[++i], "%zu", &skip) != 1 )
                usage (argv[0]);
            break;
          };
          case 'r': {
            if ( i + 2 >= argc ||
                sscanf (argv[++i], "%lf", &low) != 1 ||
                sscanf (argv[++i], "%lf", &high) != 1 )
              {
                usage (argv[0]);
              }
            break;
          }
          case 'l': {
            EXTRA_MIN = true;
            break;
          }
          case 'a': {
            if ( i + 3 >= argc ||
                sscanf (argv[++i], "%zu", &add_column) != 1 ||
                sscanf (argv[++i], "%lf", &add_exponent) != 1 ||
                sscanf (argv[++i], "%lf", &add_amplitude) != 1 )
              {
                usage (argv[0]);
              }
            break;
          };
          case 'i': {
            if ( i + 1 >= argc || sscanf (argv[++i], "%zu", &intervals) != 1 )
                usage (argv[0]);
            break;
          }
          default: {
            usage (argv[0]);
            break;
          }
        }
    }

  printf ("Range is [%f, %f], number of intervals: %zu\n", low, high, intervals);
  double incr = (high - low) / intervals;

#if defined(PREC_GMP)
  mpf_set_default_prec (256);
#endif

  ///////////////////////////////////////////////////////////////////
  // read in the data
  printf ("Opening '%s' for reading..\n", fn);
  FILE *fp = fopen (fn, "r");
  if ( !fp )
      usage (argv[0]);

  N = 0;
  M = (size_t)-1;
  while (1)
    {
      char line[4096];
      if ( !fgets (line, sizeof(line), fp) )
          break;
      if ( !isdigit (line[0]) )
          continue;

      std::istringstream stream(line);

      size_t m;
      char c;
      stream >> A[N] >> c >> m >> c >> NS[N];

      if ( m > MAXM )
        {
          std::cerr << "Data columns exceeds limit!\n";
          exit(-1);
        }

      prec_t DATA[m][m];

      for ( size_t i = 0; i < m; i++ )
        {
          stream >> c >> X[N][i];
          X[N][i] /= NS[N];
        }

      for ( size_t i = 0; i < m; i++ )
        {
          for ( size_t j = i; j < m; j++ )
            {
              stream >> c >> DATA[i][j];
              DATA[i][j] = v2d(DATA[i][j] / NS[N]);
            }
        }

      for ( size_t i = 0; i < m; i++ )
        {
          for ( size_t j = i; j < m; j++ )
            {
              COVM[N][i][j] = (DATA[i][j] - X[N][i] * X[N][j]) / NS[N];
              COVM[N][j][i] = COVM[N][i][j];
            }
        }

      M = std::min (M, m);

      if ( skip )
        {
          printf ("Skipping a record at A=%.0f\n", v2d(A[N]));
          skip--;
        }
      else
        {
          N++;
        }
    }

  fclose (fp);

  if ( N > MN )
      exit(-1);

  if ( !uv.empty() )
    {
      for ( size_t i = 1; i <= M; i++ )
        {
          if ( std::find(uv.begin(), uv.end(), i) == uv.end() )
              dv.push_back (i);
        }
    }

  std::sort (dv.begin(), dv.end(), std::greater<size_t>());
  std::stringstream dstr;
  for ( size_t i = 0; i < dv.size(); i++ )
    {
      size_t discard = dv[i];

      if ( discard > M )
        {
          printf ("Can't discard column %zu\n", discard);
          exit(-1);
        }

      dstr << " " << discard;

      for ( size_t i = 0; i < N; i++ )
        {
          for ( size_t j = discard - 1; j < M - 1; j++ )
              X[i][j] = X[i][j + 1];

          for ( size_t j = discard - 1; j < M - 1; j++ )
              for ( size_t k = 0; k < MAXM; k++ )
                  COVM[i][k][j] = COVM[i][k][j + 1];

          for ( size_t j = discard - 1; j < M - 1; j++ )
              for ( size_t k = 0; k < MAXM; k++ )
                  COVM[i][j][k] = COVM[i][j + 1][k];
        }

      for ( size_t i = 0; i < N; i++ )
        {
          X[i][M - 1] = 0;
          for ( size_t j = 0; j < M; j++ )
            {
              COVM[i][j][M - 1] = 0;
              COVM[i][M - 1][j] = 0;
            }
        }

      M--;
    }

  if ( !dstr.str().empty() )
      printf ("Discarding data columns%s\n", dstr.str().c_str());

  for ( size_t i = 0; i < fv.size(); i++ )
    {
      double x = fv[i];
      printf ("Adding a constant exponent %.4f\n", x);

      if ( M == MAXM )
        {
          fprintf (stderr, "%zu data vectors already present, more not supported!\n", M);
          exit(-1);
        }

      for ( size_t i = 0; i < N; i++ )
        {
          X[i][M] = mypow(A[i], x);
          for ( size_t j = 0; j <= M; j++ )
            {
              COVM[i][M][j] = 0;
              COVM[i][j][M] = 0;
            }
        }

      M++;
    }

  if ( add_column )
    {
      size_t k = add_column - 1;
      if ( k >= M )
        {
          fprintf (stderr, "No column %zu available!\n", add_column);
          exit(-1);
        }
      prec_t coeff = add_amplitude * sqrt(COVM[0][k][k]) / mypow(A[0], add_exponent);
      printf ("Adding %.2e*A**%.2f to all rows of column %zu (+%.4f for row 1)\n",
          v2d(coeff), v2d(add_exponent), add_column,
          v2d(coeff * mypow(A[0], add_exponent)));
      for ( size_t i = 0; i < N; i++ )
        {
          prec_t add = coeff * mypow(A[i], add_exponent);
          X[i][k] += add;
        }
    }

  if ( disp )
    {
      for ( size_t i = 0; i < N; i++ )
        {
          printf ("X:");
          for ( size_t j = 0; j < M; j++ )
              printf (" %f", v2d(X[i][j]));
          printf ("\n");
          for ( size_t j = 0; j < M; j++ )
            {
              printf ("COVM:");
              for ( size_t k = 0; k < M; k++ )
                  printf (" %.3e", v2d(COVM[i][j][k]));
              printf ("\n");
            }
        }
    }

  if ( N <= M + 1 )
    {
      fprintf (stderr, "Not enough data, N - M - 1 < 1!\n");
      exit(-1);
    }

  // calculate the inverse covariance matrices
  calc_icovm();

  size_t gc = 0;
  RM.resize(M, M);
  D.resize(M);

  ///////////////////////////////////////////////////////////////////
  // calculate the residuals and find the dimensions
  std::vector<std::pair<prec_t, prec_t> > graph;
  std::deque<std::pair<prec_t, prec_t> > last;
  for ( prec_t d = low; d <= high; d += incr )
    {
      prec_t S = calc_res (d);

      graph.push_back (std::make_pair (d, S));
      last.push_back (std::make_pair (d, S));
      if ( last.size() == 4 )
          last.pop_front();
      // find the minimums
      if ( !nomins && last.size() == 3 && last[0].second > last[1].second && last[2].second > last[1].second )
        {
          prec_t md = find_min (calc_res_functor(), last[0].first, last[1].first, last[2].first);
          prec_t y = calc_res (md);
          prec_t chi2 = chi2inv (N - M - 1);
          bool good = (y < chi2);
          prec_t search = y + chi2inv(1);

          prec_t xoff1 = find_value (search, md, incr, high) - md;
          prec_t xoff2 = md - find_value (search, md, -incr, low);
          if ( xoff1 < 0 ) xoff1 = -xoff1;
          if ( xoff2 < 0 ) xoff2 = -xoff2;
          prec_t error = std::max (xoff1, xoff2);
          printf ("[%7.4f, %7.4f]: %11.8f (y=%.3e), err=%.8f, g=%d\n",
              v2d(last[0].first), v2d(last[2].first), v2d(md), v2d (y),
              v2d(error), good ? 1 : 0);
          showC(md);
          if ( apply_traditional && gc < M )
            {
              calc_res(md, &RM(gc, 0));
              D[gc] = md;
              gc++;
            }

          graph.push_back (std::make_pair (md, y));
        }
    }

  if ( apply_traditional )
    {
      RM.invert();
#if 1
    if ( gc != M ) {
      for ( size_t i = 0; i < M; i++ )
          D[i] = i;
      //D[0] = 0;
      /*
      D[0] = -1;
      D[1] = 0;
      D[2] = 1.39;
      */
      for ( size_t r = 0; r < M; r++ )
          for ( size_t c = 0; c < M; c++ )
              RM(r, c) = c == 0 ? -1.0 : 1.0;
      }
#endif
      RM.print();

      prec_t chi2level = chi2inv(N*M-M*M-M);
      prec_t chi2 = calc_chi2overall();
      printf ("Overall chi2=%g (good=%d,%g), df=%zu\n", v2d(chi2), chi2 < chi2level ? 1 : 0, v2d(chi2level), N*M-M*M-M);

      make_alpha_beta();

      prec_t old_chi;
      std::vector<prec_t> old_D;
      matrix old_RM;
      prec_t lambda = 1e-6;
      matrix da;

      while (lambda < 1e9)
        {
          old_chi = calc_chi2overall();
          old_D = D;
          old_RM = RM;
          make_alpha_beta();
          for ( size_t r = 0; r < M; r++ )
              alpha(r, r) *= (1 + lambda);
          alpha.invert();
          matrix::mul(da, alpha, beta);
          size_t dai = 0;
          for ( size_t i = 0; i < M; i++ )
            {
              D[i] += da(dai, 0);
              dai++;
            }
          for ( size_t r = 0; r < M; r++ )
            {
              for ( size_t c = 0; c < M; c++ )
                {
                  RM(r,c) += da(dai, 0);
                  dai++;
                }
            }
          prec_t new_chi = calc_chi2overall();
          if ( old_chi < new_chi )
            {
              RM.swap(old_RM);
              D.swap(old_D);
              lambda = lambda * 10;
              continue;
            }
          if ( lambda > 1e-16 )
              lambda = lambda / 10;
          printf("Loop %g, D=", v2d(calc_chi2overall()));
          for ( size_t i = 0; i < M; i++ )
              printf (" %g", v2d(D[i]));

          make_alpha_beta();
          alpha.invert();
          printf(" +- %g\n", v2d(1.96 * sqrt(alpha(M-1, M-1))));

          if ( myabs(old_chi - new_chi) < 1e-8 )
              break;
        }

      RM.print();

      chi2 = calc_chi2overall();
      printf ("Overall chi2=%g (good=%d,%g), df=%zu\n", v2d(chi2), chi2 < chi2level ? 1 : 0, v2d(chi2level), N*M-M*M-M);
      make_alpha_beta();
      alpha.invert();
      for ( size_t i = 0; i < M; i++ )
          printf ("D%zu =  %g +- %g\n", i, v2d(D[i]), v2d(1.96 * sqrt(alpha(i, i))));
    }

#if 0
  {
    FILE *tfp = fopen("graph.csv", "w");
    for ( size_t i = 0; i < N; i++ )
      {
        fprintf (tfp, "%g", v2d(A[i]));
        for ( size_t j = 0; j < M; j++ )
          {
            prec_t res = 0;
            for ( size_t r = 0; r < M; r++ )
                res += RM(j, r)*mypow(A[i], D[r]);
            fprintf (tfp, ",%g,%g", v2d(X[i][j]), v2d(res));
          }
        fprintf (tfp, "\n");
      }
    fclose (tfp);
  }
#endif

  printf ("Writing graph data into 'calc.out'\n");
  FILE *outfp = fopen ("calc.out", "w");
  if ( !outfp )
      exit(-1);
  std::sort (graph.begin(), graph.end());
  for ( size_t i = 0; i < graph.size(); i++ )
      fprintf (outfp, "%.16f %.16e\n", v2d(graph[i].first), v2d(graph[i].second));
  fclose (outfp);
}

